Introduction to Bayesian statistical modelling

A course with R, Stan, and brms

Ladislas Nalborczyk (UNICOG, NeuroSpin, CEA, Gif/Yvette, France)

Planning

Course n°01: Introduction to Bayesian inference, Beta-Binomial model
Course n°02: Introduction to brms, linear regression
Course n°03: Markov Chain Monte Carlo, generalised linear model
Course n°04: Multilevel models, cognitive models

\[\newcommand\given[1][]{\:#1\vert\:}\]

Multilevel models

The aim is to build a model that can learn at several levels, that can produce estimates that will be informed by the different groups present in the data. We will follow the following example throughout this course.

Let’s imagine that we’ve built a robot that visits cafés, and that has fun measuring the waiting time after ordering. This robot visits 20 different cafés, 5 times in the morning and 5 times in the afternoon, and measures the time (in minutes) it takes to get a coffee.

Coffee robot

library(tidyverse)
library(imsb)

df <- open_data(robot)
head(x = df, n = 15)
   cafe afternoon      wait
1     1         0 4.9989926
2     1         1 2.2133944
3     1         0 4.1866730
4     1         1 3.5624399
5     1         0 3.9956779
6     1         1 2.8957176
7     1         0 3.7804582
8     1         1 2.3844837
9     1         0 3.8617982
10    1         1 2.5800004
11    2         0 2.7421223
12    2         1 1.3525907
13    2         0 2.5215095
14    2         1 0.9628102
15    2         0 1.9543977

Coffee robot

df %>%
  ggplot(aes(x = factor(cafe), y = wait, fill = factor(afternoon) ) ) +
  geom_dotplot(
    stackdir = "center", binaxis = "y",
    dotsize = 1, show.legend = FALSE
    ) +
  geom_hline(yintercept = mean(df$wait), linetype = 3) +
  facet_wrap(~afternoon, ncol = 2) +
  labs(x = "Café", y = "Waiting time (min)")

Coffee robot, a first model

An initial model can be built, estimating the average time (across all cafés combined) to be served.

\[ \begin{align} \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(5, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{HalfCauchy}(0, 2)} \\ \end{align} \]

Half-Cauchy

\[ p(x \given x_{0}, \gamma) = \left(\pi \gamma \left[1 + \left(\frac{x-x_{0}}{\gamma}\right)^{2}\right] \right)^{-1} \]

ggplot(data = data.frame(x = c(0, 10) ), aes(x = x) ) +
    stat_function(
        fun = dcauchy,
        args = list(location = 0, scale = 2), size = 1.5
        )

Coffee robot, a first model

library(brms)

mod1 <- brm(
  formula = wait ~ 1,
  prior = c(
    prior(normal(5, 10), class = Intercept),
    prior(cauchy(0, 2), class = sigma)
    ),
  data = df,
  cores = parallel::detectCores()
  )
posterior_summary(x = mod1, probs = c(0.025, 0.975), pars = c("^b_", "sigma") )
            Estimate  Est.Error     Q2.5    Q97.5
b_Intercept 3.121625 0.07923713 2.966771 3.280714
sigma       1.143018 0.05640463 1.036736 1.256459

Diagnostic plot

plot(
  x = mod1, combo = c("dens_overlay", "trace"),
  theme = theme_bw(base_size = 16, base_family = "Open Sans")
  )

One intercept per café

Second model which estimates one intercept per coffee. Equivalent to constructing 20 dummy variables.

\[ \begin{align} \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha_{\text{café}[i]}} \\ \color{steelblue}{\alpha_{\text{café}[i]}} \ &\color{steelblue}{\sim \mathrm{Normal}(5, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{HalfCauchy}(0, 2)} \\ \end{align} \]

mod2 <- brm(
  formula = wait ~ 0 + factor(cafe),
  prior = c(
    prior(normal(5, 10), class = b),
    prior(cauchy(0, 2), class = sigma)
    ),
  data = df,
  cores = parallel::detectCores()
  )

One intercept per café

posterior_summary(x = mod2, pars = "^b_")
               Estimate Est.Error      Q2.5    Q97.5
b_factorcafe1  3.446599 0.2674826 2.9112647 3.967342
b_factorcafe2  1.730125 0.2615548 1.2190581 2.238199
b_factorcafe3  3.321925 0.2629713 2.8148474 3.833974
b_factorcafe4  2.803049 0.2584659 2.2855814 3.315254
b_factorcafe5  1.464544 0.2627843 0.9555137 1.969844
b_factorcafe6  3.643927 0.2711427 3.1090742 4.194942
b_factorcafe7  2.944457 0.2614577 2.4278226 3.472341
b_factorcafe8  3.172552 0.2660662 2.6328560 3.694486
b_factorcafe9  3.337923 0.2689664 2.8118534 3.865261
b_factorcafe10 3.099324 0.2625124 2.5673710 3.619401
b_factorcafe11 1.916332 0.2571375 1.4105794 2.419089
b_factorcafe12 3.494184 0.2541949 3.0054751 3.998191
b_factorcafe13 3.220359 0.2545903 2.7332581 3.714874
b_factorcafe14 2.626250 0.2633425 2.1031520 3.147996
b_factorcafe15 3.476897 0.2641088 2.9619580 3.991096
b_factorcafe16 2.997862 0.2599787 2.4853711 3.510518
b_factorcafe17 3.884514 0.2633627 3.3741596 4.398071
b_factorcafe18 5.529265 0.2612210 5.0195458 6.061847
b_factorcafe19 2.977062 0.2605880 2.4784251 3.489196
b_factorcafe20 3.364699 0.2530715 2.8744273 3.870453

Multilevel model

Couldn’t we ensure that the time measured at café 1 informs the measurement taken at coffee 2 and coffee 3? As well as the average time taken to be served? We’re going to learn the priors from the data…

\[ \begin{align} \text{Level 1}: \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha_{\text{café}[i]}} \\ \text{Level 2}: \color{steelblue}{\alpha_{\text{café}}} \ &\color{steelblue}{\sim \mathrm{Normal}(\alpha,\sigma_{\text{café}})} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(5, 10)} \\ \color{steelblue}{\sigma_{\text{café}}} \ &\color{steelblue}{\sim \mathrm{HalfCauchy}(0, 2)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{HalfCauchy}(0, 2)} \\ \end{align} \]

The prior for the intercept of each coffee (\(\alpha_{\text{café}}\)) is now a function of two parameters (\(\alpha\) and \(\sigma_{\text{café}}\)). \(\alpha\) and \(\sigma_{\text{café}}\) are called hyper-parameters, they are parameters for parameters, and their priors are called hyperpriors. There are two levels in the model…

Equivalences

\[ \begin{align} \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha_{\text{café}[i]}} \\ \color{steelblue}{\alpha_{\text{café}}} \ &\color{steelblue}{\sim \mathrm{Normal}(\alpha,\sigma_{\text{café}})} \\ \end{align} \]

NB: \(\alpha\) is defined here in the prior for \(\alpha_{\text{café}}\) but it could, in the same way, be defined in the linear model:

\[ \begin{align} \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \alpha_{\text{café}[i]}} \\ \color{steelblue}{\alpha_{\text{café}}} \ &\color{steelblue}{\sim \mathrm{Normal}(0,\sigma_{\text{café}})} \\ \end{align} \]

We can always “remove” the mean from a Gaussian distribution and consider it as a constant plus a Gaussian centred on zero.

NB: when \(\alpha\) is defined in the linear model, the \(\alpha_{\text{café}}\) represent deviations from the mean intercept. It is therefore necessary to add \(\alpha\) and \(\alpha_{\text{café}}\) to obtain the average waiting time per café…

Equivalences

y1 <- rnorm(n = 1e4, mean = 5, sd = 1)
y2 <- rnorm(n = 1e4, mean = 0, sd = 1) + 5

data.frame(y1 = y1, y2 = y2) %>%
    pivot_longer(cols = 1:2, names_to = "x", values_to = "y") %>%
    ggplot(aes(x = y, colour = x) ) +
    geom_density(show.legend = FALSE)

Multilevel model

mod3 <- brm(
  formula = wait ~ 1 + (1 | cafe),
  prior = c(
    prior(normal(5, 10), class = Intercept),
    prior(cauchy(0, 2), class = sigma),
    prior(cauchy(0, 2), class = sd)
    ),
  data = df,
  warmup = 1000, iter = 5000,
  cores = parallel::detectCores()
  )

This model has 23 parameters, the general intercept \(\alpha\), the residual variability \(\sigma\), the variability between cafés and one intercept per café.

Shrinkage

Shrinkage magic (Efron & Morris, 1977)

L’estimateur James-Stein est défini comme \(z = \bar{y} + c(y - \bar{y})\), où \(\bar{y}\) désigne la moyenne de l’échantillon, \(y\) une observation individuelle, et \(c\) une constante, le shrinking factor (Efron & Morris, 1977).

Shrinkage magic (Efron & Morris, 1977)

Le shrinking factor est déterminé à la fois par la variabilité (imprécision) de la mesure (e.g., son écart-type) et par la distance à l’estimation moyenne (i.e., \(y - \bar{y}\)). En d’autres termes, cet estimateur fait moins “confiance” (i.e., accorde moins de poids) aux observations imprécises et/ou extrêmes. En pratique, le shrinkage agit comme une protection contre le sur-apprentissage (overfitting).

Pooling

Le shrinkage observé slide précédente est dû à des phénomènes de partage (pooling) de l’information entre les cafés. L’estimation de l’intercept pour chaque café informe l’estimation de l’intercept des autres cafés, ainsi que l’estimation de l’intercept général (i.e., la moyenne générale).

On distingue en général trois perspectives (ou stratégies) :

  • Complete pooling : on suppose que le temps d’attente est invariant, on estime un intercept commun (mod1).

  • No pooling : on suppose que les temps d’attente de chaque café sont uniques et indépendants : on estime un intercept par café, mais sans informer le niveau supérieur (mod2).

  • Partial pooling : on utilise un prior adaptatif, comme dans l’exemple précédent (mod3).

La stratégie complete pooling en général underfitte les données (faibles capacités de prédiction) tandis que le stratégie no pooling revient à overfitter les données (faibles capacités de prédiction ici aussi). La stratégie partial pooling (i.e., celle des modèles multi-niveaux) permet d’’équilibrer underfitting et overfitting.

Comparaison de modèles

On peut comparer ces trois modèles en utilisant le WAIC (discuté au Cours n°07).

# calcul du WAIC et ajout du WAIC à chaque modèle
mod1 <- add_criterion(mod1, "waic")
mod2 <- add_criterion(mod2, "waic")
mod3 <- add_criterion(mod3, "waic")

# comparaison des WAIC de chaque modèle
w <- loo_compare(mod1, mod2, mod3, criterion = "waic")
print(w, simplify = FALSE)
     elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic   se_waic
mod3    0.0       0.0  -253.9       8.3         18.3    1.5     507.8   16.6 
mod2   -0.9       1.3  -254.8       8.4         19.8    1.6     509.6   16.8 
mod1  -57.2      10.6  -311.1      10.5          2.0    0.3     622.1   21.1 

On remarque que le modèle 3 a seulement 18 “effective parameters” (pWAIC) et moins de paramètres que le modèle 2, alors qu’il en a en réalité 2 de plus… posterior_summary(mod3)[3, 1] nous donne le sigma du prior adaptatif des \(\alpha_{\text{café}}\) (\(\sigma_{\text{café}} = 0.82\)). On remarque que ce sigma est très faible et correspond à assigner un prior très contraignant, ou régularisateur.

Comparaison de modèles

On compare les estimations du premier modèle (complete pooling model) et du troisième modèle (partial pooling model).

posterior_summary(mod1, pars = c("^b", "sigma") )
            Estimate  Est.Error     Q2.5    Q97.5
b_Intercept 3.121625 0.07923713 2.966771 3.280714
sigma       1.143018 0.05640463 1.036736 1.256459
posterior_summary(mod3, pars = c("^b", "sigma") )
             Estimate  Est.Error     Q2.5     Q97.5
b_Intercept 3.1247949 0.20583672 2.713431 3.5217680
sigma       0.8215185 0.04361643 0.741872 0.9129304

Les deux modèles font la même prédiction (en moyenne) pour \(\alpha\), mais le modèle 3 est plus incertain de sa prédiction que le modèle 1 (voir l’erreur standard pour \(\alpha\))…

L’estimation de \(\sigma\) du modèle 3 est plus petite que celle du modèle 1 car le modèle 3 décompose la variabilité non expliquée en deux sources : la variabilité du temps d’attente entre les cafés \(\sigma_{\text{café}}\) et la variabilité résiduelle \(\sigma\).

Robot et café

Imaginons que notre robot ne visite pas tous les cafés le même nombre de fois (comme dans le cas précédent) mais qu’il visite plus souvent les cafés proches de chez lui…

df2 <- open_data(robot_unequal) # nouveau jeu de données

mod4 <- brm(
  formula = wait ~ 1 + (1 | cafe),
  prior = c(
    prior(normal(5, 10), class = Intercept),
    prior(cauchy(0, 2), class = sigma),
    prior(cauchy(0, 2), class = sd)
    ),
  data = df2,
  warmup = 1000, iter = 5000,
  cores = parallel::detectCores()
  )

Shrinkage

On observe que les cafés qui sont souvent visités (à droite) subissent moins l’effet du shrinkage. Leur estimation est moins “tirée” vers la moyenne générale que les estimations des cafés les moins souvent visités (à gauche).

Aparté : effets fixes et effets aléatoires

Cinq définitions (contradictoires) relevées par Gelman (2005).

  • Fixed effects are constant across individuals, and random effects vary.
  • Effects are fixed if they are interesting in themselves or random if there is interest in the underlying population.
  • When a sample exhausts the population, the corresponding variable is fixed; when the sample is a small (i.e., negligible) part of the population the corresponding variable is random.
  • If an effect is assumed to be a realized value of a random variable, it is called a random effect.
  • Fixed effects are estimated using least squares (or, more generally, maximum likelihood) and random effects are estimated with shrinkage.

Gelman & Hill (2006) suggèrent plutôt l’utilisation des termes de constant effects et varying effects, et de toujours utiliser la modélisation multi-niveaux, en considérant que ce qu’on appelle effet fixe peut simplement être considéré comme un effet aléatoire dont la variance serait égale à \(0\).

Régularisation et terminologie

Le fait de faire varier les intercepts de chaque café est simplement une autre manière de régulariser (de manière adaptative), c’est à dire de diminuer le poids accordé aux données dans l’estimation. Le modèle devient à même d’estimer à quel point les groupes (ici les cafés) sont différents, tout en estimant les caractéristiques de chaque café…

Différence entre les cross-classified (ou “crossed”) multilevel models et nested or hierarchical multilevel models. Le premier type de modèle concerne des données structurées selon deux (ou plus) facteurs aléatoires non “nichés”. Le deuxième type de modèles concerne des données structurées de manière hiérarchique (e.g., un élève dans une classe dans une école dans une ville…). Voir cette discussion pour plus de détails.

Les deux types de modèles s’écrivent cependant de manière similaire, sur plusieurs “niveaux”. Le terme “multi-niveaux” (dans notre terminologie) fait donc référence à la structure du modèle, à sa spécification. À distinguer de la structure des données.

Exemple de modèle “cross-classified”

On pourrait se poser la question de savoir si la récence des cafés (leur âge) ne serait pas une source de variabilité non contrôlée ? Il suffit d’ajouter un intercept qui varie par âge, et de lui attribuer un prior adaptatif.

\[ \begin{align} \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \alpha_{\text{café}[i]} + \alpha_{\text{âge}[i]}} \\ \color{steelblue}{\alpha_{\text{café}}} \ &\color{steelblue}{\sim \mathrm{Normal}(5, \sigma_{\text{café}})} \\ \color{steelblue}{\alpha_{\text{âge}}} \ &\color{steelblue}{\sim \mathrm{Normal}(5, \sigma_{\text{âge}})} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\sigma_{\text{café}}} \ &\color{steelblue}{\sim \mathrm{HalfCauchy}(0, 2)} \\ \color{steelblue}{\sigma_{\text{âge}}} \ &\color{steelblue}{\sim \mathrm{HalfCauchy}(0, 2)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{HalfCauchy}(0, 2)} \\ \end{align} \]

Robot et café : varying intercept + varying slope

On s’intéresse maintenant à l’effet du moment de la journée sur le temps d’attente. Attend-on plus le matin, ou l’après-midi ?

\[ \begin{align} \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha_{\text{café}[i]} + \beta_{\text{café}[i]} \times A_{i}} \\ \end{align} \]

\(A_{i}\) est une dummy variable codée 0/1 pour le matin et l’après-midi et où \(\beta_{\text{café}}\) est donc un paramètre de différence (i.e., une pente) entre le matin et l’après-midi.

Remarque : on sait que les cafés ont des intercepts et des pentes qui co-varient… Les cafés populaires seront surchargés le matin et beaucoup moins l’après-midi, résultant en une pente importante. Ces cafés auront aussi un temps d’attente moyen plus long (i.e., un intercept plus grand). Dans ces cafés, \(\alpha\) est grand et \(\beta\) est loin de zéro. À l’inverse, dans un café peu populaire, le temps d’attente sera faible, ainsi que la différence entre matin et après-midi.

On pourrait donc utiliser la co-variation entre intercept et pente pour faire de meilleures inférences. Autrement dit, faire en sorte que l’estimation de l’intercept informe celle de la pente, et réciproquement.

Robot et café : varying intercept + varying slope

On s’intéresse maintenant à l’effet du moment de la journée sur le temps d’attente. Attend-on plus le matin, ou l’après-midi ?

\[ \begin{align} \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha_{\text{café}[i]} + \beta_{\text{café}[i]} \times A_{i}} \\ \color{steelblue}{\begin{bmatrix} \alpha_{\text{café}} \\ \beta_{\text{café}} \\ \end{bmatrix}} \ &\color{steelblue}{\sim \mathrm{MVNormal}\bigg(\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \textbf{S}\bigg)} \\ \end{align} \]

La troisième ligne postule que chaque café a un intercept \(\alpha_{\text{café}}\) et une pente \(\beta_{\text{café}}\), définis par un prior Gaussien bivarié (i.e., à deux dimensions) ayant comme moyennes \(\alpha\) et \(\beta\) et comme matrice de covariance \(\textbf{S}\).

Aparté : distribution gaussienne multivariée

\[\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})\]

\(\boldsymbol{\mu}\) est un vecteur (à \(k\) dimensions) de moyennes, par exemple: mu <- c(a, b).

\(\boldsymbol{\Sigma}\) est une matrice de covariance de \(k \times k\) dimensions, et qui correspond à la matrice donnée par la fonction vcov().

\[ \begin{align} \boldsymbol{\Sigma} &= \begin{pmatrix} \sigma_{\alpha}^{2} & \sigma_{\alpha} \sigma_{\beta} \rho \\ \sigma_{\alpha} \sigma_{\beta} \rho & \sigma_{\beta}^{2} \\ \end{pmatrix} \\ \end{align} \]

Aparté : distribution gaussienne multivariée

Aparté : distribution gaussienne multivariée

\[ \begin{align} \boldsymbol{\Sigma} &= \begin{pmatrix} \sigma_{\alpha}^{2} & \sigma_{\alpha} \sigma_{\beta} \rho \\ \sigma_{\alpha} \sigma_{\beta} \rho & \sigma_{\beta}^{2} \\ \end{pmatrix} \\ \end{align} \]

Cette matrice peut se construire de deux manières différentes, strictement équivalentes.

sigma_a <- 1
sigma_b <- 0.75
rho <- 0.7
cov_ab <- sigma_a * sigma_b * rho
(Sigma1 <- matrix(c(sigma_a^2, cov_ab, cov_ab, sigma_b^2), ncol = 2) )
      [,1]   [,2]
[1,] 1.000 0.5250
[2,] 0.525 0.5625

Aparté : distribution gaussienne multivariée

\[ \begin{align} \boldsymbol{\Sigma} &= \begin{pmatrix} \sigma_{\alpha}^{2} & \sigma_{\alpha} \sigma_{\beta} \rho \\ \sigma_{\alpha} \sigma_{\beta} \rho & \sigma_{\beta}^{2} \\ \end{pmatrix} \\ \end{align} \]

La deuxième méthode est pratique car elle considère séparément les écart-types et les corrélations.

(sigmas <- c(sigma_a, sigma_b) ) # standard deviations
[1] 1.00 0.75
(Rho <- matrix(c(1, rho, rho, 1), nrow = 2) ) # correlation matrix
     [,1] [,2]
[1,]  1.0  0.7
[2,]  0.7  1.0
(Sigma2 <- diag(sigmas) %*% Rho %*% diag(sigmas) )
      [,1]   [,2]
[1,] 1.000 0.5250
[2,] 0.525 0.5625

Robot et café : varying intercept + varying slope

\[ \begin{align} \color{orangered}{w_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha_{\text{café}[i]} + \beta_{\text{café}[i]} \times A_{i}} \\ \color{steelblue}{\begin{bmatrix} \alpha_{\text{café}} \\ \beta_{\text{café}} \\ \end{bmatrix}} \ &\color{steelblue}{\sim \mathrm{MVNormal}\bigg(\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \textbf{S}\bigg)} \\ \color{black}{\textbf{S}} \ &\color{black}{= \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \\ \end{pmatrix} \ \textbf{R} \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \\ \end{pmatrix}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal} (0, 10)} \\ \color{steelblue}{\beta} \ &\color{steelblue}{\sim \mathrm{Normal} (0, 10)} \\ \color{steelblue}{\sigma_{\alpha}} \ &\color{steelblue}{\sim \mathrm{HalfCauchy} (0, 2)} \\ \color{steelblue}{\sigma_{\beta}} \ &\color{steelblue}{\sim \mathrm{HalfCauchy} (0, 2)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{HalfCauchy} (0, 2)} \\ \color{steelblue}{\textbf{R}} \ &\color{steelblue}{\sim \mathrm{LKJ}(2)} \\ \end{align} \]

\(\textbf{S}\) est définie en factorisant \(\sigma_{\alpha}\), \(\sigma_{\beta}\), et la matrice de corrélation \(\textbf{R}\). La suite du modèle définit simplement les priors pour les effets constants. La dernière ligne spécifie le prior pour \(\textbf{R}\).

LKJ prior

Prior proposé par Lewandowski et al. (2009). Un seul paramètre \(\zeta\) (zeta) spécifie la concentration de la distribution du coefficient de corrélation. Le prior \(\mathrm{LKJ}(2)\) définit un prior peu informatif pour \(\rho\) (rho) qui est sceptique des corrélations extrêmes (i.e., des valeurs proches de \(-1\) ou \(1\)).

Rappels de syntaxe

Le paquet brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le paquet lme4.

Reaction ~ Days + (1 + Days | Subject)

La partie gauche représente notre variable dépendante (ou “outcome”, i.e., ce qu’on essaye de prédire).

La partie droite permet de définir les prédicteurs. L’intercept est généralement implicite, de sorte que les deux écritures ci-dessous sont équivalentes.

Reaction ~ Days + (1 + Days | Subject)
Reaction ~ 1 + Days + (1 + Days | Subject)

Rappels de syntaxe

La première partie de la partie droite de la formule représente les effets constants (effets fixes), tandis que la seconde partie (entre parenthèses) représente les effets “variants” ou “variables” (effets aléatoires).

Reaction ~ 1 + Days + (1 | Subject)
Reaction ~ 1 + Days + (1 + Days | Subject)

Le premier modèle ci-dessus contient seulement un intercept variable, qui varie par Subject. Le deuxième modèle contient également un intercept variable, mais aussi une pente variable pour l’effet de Days.

Rappels de syntaxe

Lorsqu’on inclut plusieurs effets variants (e.g., un intercept et une pente variables), brms postule qu’on souhaite aussi estimer la corrélation entre ces deux effets. Dans le cas contraire, on peut supprimer cette corrélation (i.e., la fixer à 0) en utilisant ||.

Reaction ~ Days + (1 + Days || Subject)

Les modèles précédents postulaient un modèle génératif Gaussien. Ce postulat peut être changé facilement en spécifiant la fonction souhaitée via l’argument family.

brm(formula = Reaction ~ 1 + Days + (1 + Days | Subject), family = lognormal() )

Modèle brms

On spécifie un intercept et une pente (pour l’effet d’afternoon) qui varient par cafe.

mod5 <- brm(
  formula = wait ~ 1 + afternoon + (1 + afternoon | cafe),
  prior = c(
    prior(normal(0, 10), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 2), class = sigma),
    prior(cauchy(0, 2), class = sd)
    ),
  data = df,
  warmup = 1000, iter = 5000,
  cores = parallel::detectCores()
  )

Distribution postérieure

post <- as_draws_df(x = mod5) # extracts posterior samples
R <- rethinking::rlkjcorr(n = 16000, K = 2, eta = 2) # samples from prior

data.frame(prior = R[, 1, 2], posterior = post$cor_cafe__Intercept__afternoon) %>%
    gather(type, value, prior:posterior) %>%
    ggplot(aes(x = value, color = type, fill = type) ) +
    geom_histogram(position = "identity", alpha = 0.2) +
    labs(x = expression(rho), y = "Nombre d'échantillons")

Shrinkage en deux dimensions

Comparaison de modèles

On compare le premier modèle (complete pooling model), le troisième modèle (partial pooling model), et le dernier modèle (avec intercept et pente variable).

# comparaison des WAIC de chaque modèle
mod5 <- add_criterion(mod5, "waic")
w <- loo_compare(mod1, mod2, mod3, mod5, criterion = "waic")
print(w, simplify = FALSE)
     elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic   se_waic
mod5    0.0       0.0  -155.1      10.0         26.4    2.6     310.3   20.1 
mod3  -98.8       8.3  -253.9       8.3         18.3    1.5     507.8   16.6 
mod2  -99.7       8.3  -254.8       8.4         19.8    1.6     509.6   16.8 
mod1 -155.9      13.7  -311.1      10.5          2.0    0.3     622.1   21.1 
model_weights(mod1, mod2, mod3, mod5, weights = "waic")
        mod1         mod2         mod3         mod5 
1.880249e-68 5.080298e-44 1.273836e-43 1.000000e+00 

Comparaison de modèles

L’estimation du temps d’attente moyen est plus incertaine lorsqu’on prend en compte de nouvelles sources d’erreur. Cependant, l’erreur du modèle (i.e., ce qui n’est pas expliqué), la variation résiduelle \(\sigma\), diminue…

posterior_summary(mod1, pars = c("^b", "sigma") )
            Estimate  Est.Error     Q2.5    Q97.5
b_Intercept 3.121625 0.07923713 2.966771 3.280714
sigma       1.143018 0.05640463 1.036736 1.256459
posterior_summary(mod3, pars = c("^b", "sigma") )
             Estimate  Est.Error     Q2.5     Q97.5
b_Intercept 3.1247949 0.20583672 2.713431 3.5217680
sigma       0.8215185 0.04361643 0.741872 0.9129304
posterior_summary(mod5, pars = c("^b", "sigma") )
              Estimate  Est.Error       Q2.5      Q97.5
b_Intercept  3.7355470 0.21026789  3.3160121  4.1495370
b_afternoon -1.2319627 0.08702830 -1.4031316 -1.0608732
sigma        0.4898022 0.02720104  0.4394355  0.5462991

Conclusions

Les modèles multi-niveaux (ou “modèles mixtes”) sont des extensions naturelles des modèles de régression classiques, où les paramètres de ces derniers se voient eux-même attribués des “modèles”, gouvernés par des hyper-paramètres.

Cette extension permet de faire des prédictions plus précises en prenant en compte la variabilité liée aux groupes ou structures (clusters) présent(e)s dans les données. Autrement dit, en modélisant les populations d’où sont tirés les effets aléatoires (e.g., la population de participants ou de stimuli).

Un modèle de régression classique est équivalent à un modèle multi-niveaux où la variabilité des effets aléatoires serait fixée à \(0\).

La cadre bayésien permet une interprétation naturelle des distributions desquelles proviennent les effets aléatoires (varying effects). En effet, ces distributions peuvent être interprétées comme des distributions a priori, dont les paramètres sont estimés à partir des données.

Travaux pratiques - sleepstudy

library(lme4)
data(sleepstudy)
head(sleepstudy, 20)
   Reaction Days Subject
1  249.5600    0     308
2  258.7047    1     308
3  250.8006    2     308
4  321.4398    3     308
5  356.8519    4     308
6  414.6901    5     308
7  382.2038    6     308
8  290.1486    7     308
9  430.5853    8     308
10 466.3535    9     308
11 222.7339    0     309
12 205.2658    1     309
13 202.9778    2     309
14 204.7070    3     309
15 207.7161    4     309
16 215.9618    5     309
17 213.6303    6     309
18 217.7272    7     309
19 224.2957    8     309
20 237.3142    9     309

Travaux pratiques - sleepstudy

sleepstudy %>%
    ggplot(aes(x = Days, y = Reaction) ) +
    geom_smooth(method = "lm", colour = "black") +
    geom_point() +
    facet_wrap(~Subject, nrow = 2) +
    scale_x_continuous(breaks = c(0, 2, 4, 6, 8) )

Travaux pratiques - sleepstudy

À vous de construire les modèles mathématiques et les modèles brms correspondant aux modèles suivants :

  • Modèle avec seulement l’effet fixe de Days.
  • Modèle avec l’effet fixe de Days + un effet aléatoire de Subject (varying intercept).
  • Modèle avec l’effet fixe de Days + un effet aléatoire de Subject. (varying intercept) + un effet aléatoire de Days (varying slope).

Comparez ensuite ces modèles en utilisant les outils discutés aux cours précédents (e.g., WAIC) et concluez.

Proposition de solution

fmod0 <- lm(Reaction ~ Days, sleepstudy)
fmod1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
fmod2 <- lmer(Reaction ~ Days + (1 + Days | Subject), sleepstudy)

anova(fmod1, fmod2)
Data: sleepstudy
Models:
fmod1: Reaction ~ Days + (1 | Subject)
fmod2: Reaction ~ Days + (1 + Days | Subject)
      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
fmod1    4 1802.1 1814.8 -897.04   1794.1                         
fmod2    6 1763.9 1783.1 -875.97   1751.9 42.139  2  7.072e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Proposition de solution

mod6 <- brm(
  Reaction ~ 1 + Days,
  prior = c(
    prior(normal(200, 100), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sigma)
    ),
  data = sleepstudy,
  warmup = 1000, iter = 5000,
  cores = parallel::detectCores()
  )
posterior_summary(mod6)
              Estimate Est.Error        Q2.5      Q97.5
b_Intercept  251.97322 6.5404674  239.347696  264.75200
b_Days        10.31388 1.2345721    7.930818   12.72380
sigma         47.80875 2.5655356   43.137623   53.22188
lprior       -15.69346 0.1653433  -16.043111  -15.38961
lp__        -963.47490 1.2329454 -966.725660 -962.07522

Proposition de solution

mod7 <- brm(
  Reaction ~ 1 + Days + (1 | Subject),
  prior = c(
    prior(normal(200, 100), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sigma),
    prior(cauchy(0, 10), class = sd)
    ),
  data = sleepstudy,
  warmup = 1000, iter = 5000,
  cores = parallel::detectCores()
  )
posterior_summary(mod7, pars = c("^b", "sigma") )
             Estimate Est.Error       Q2.5    Q97.5
b_Intercept 250.89457 9.9670291 231.292797 270.6747
b_Days       10.38816 0.7907073   8.839429  11.9414
sigma        31.05272 1.7448405  27.876041  34.6853

Proposition de solution

mod8 <- brm(
  Reaction ~ 1 + Days + (1 + Days | Subject),
  prior = c(
    prior(normal(200, 100), class = Intercept),
    prior(normal(0, 10), class = b),
    prior(cauchy(0, 10), class = sigma),
    prior(cauchy(0, 10), class = sd)
    ),
  data = sleepstudy,
  warmup = 1000, iter = 5000,
  cores = parallel::detectCores()
  )
posterior_summary(mod8, pars = c("^b", "sigma") )
             Estimate Est.Error       Q2.5     Q97.5
b_Intercept 250.96896  6.949501 237.336539 264.81910
b_Days       10.10206  1.626173   6.893363  13.26584
sigma        25.88374  1.586442  23.006100  29.26188

Proposition de solution

# calcul du WAIC et ajout du WAIC à chaque modèle
mod6 <- add_criterion(mod6, "waic")
mod7 <- add_criterion(mod7, "waic")
mod8 <- add_criterion(mod8, "waic")

# comparaison des WAIC de chaque modèle
w <- loo_compare(mod6, mod7, mod8, criterion = "waic")
print(w, simplify = FALSE)
     elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic   se_waic
mod8    0.0       0.0  -860.5      22.3         33.0    8.4    1721.0   44.6 
mod7  -24.2      11.5  -884.7      14.5         19.4    3.3    1769.5   28.9 
mod6  -92.8      20.9  -953.3      10.6          3.2    0.5    1906.6   21.1 
# calcul du poids relatif de chaque modèle
model_weights(mod6, mod7, mod8, weights = "waic")
        mod6         mod7         mod8 
5.015604e-41 2.977648e-11 1.000000e+00 

Bayesian workflow (Gelman et al., 2020)

Scientific and cognitive modelling

What’s a model?

What’s in a model?

Assumptions (not arbitrary!)…

What is a model good for?

“One of the most basic problem in scientific inference is the so-called inverse problem: How to figure out causes from observations. It is a problem, because many different causes can produce the same evidence. So while it can be easy to go forward from a known cause to predicted observations, it can very hard to go backwards from observation to cause” (McElreath, 2020)…

Geometric people

Evidence accumulation models

Drift diffusion model

Conclusions

La statistique bayésienne est une approche générale de l’estimation de paramètres. Cette approche utilise la théorie des probabilités pour quantifier l’incertitude vis à vis de la valeur des paramètres de modèles statistiques.

Ces modèles sont composés de différents blocs (e.g., fonction de vraisemblance, priors, modèle linéaire ou non-linéaire) qui sont modifiables à souhait. Ce qu’on appelle classiquement “conditions d’application” sont simplement les conséquences des choix de modélisation réalisés par l’utilisateur. Autrement dit, c’est l’utilisateur qui choisit (et ne subit pas) les conditions d’application.

Nous avons vu que le modèle de régression linéaire est un modèle très flexible qui permet de décrire, via la modification de la fonction de vraisemblance et via l’introduction de fonctions de lien, des relations complexes (e.g., non-linéaires) entre variable prédite et variables prédictrices. Ces modèles peuvent gagner en précision par la prise en compte de la variabilité et des structures présentes dans les données (cf. modèles multi-niveaux).

Conclusions

Le paquet brms est un véritable couteau suisse de l’analyse statistique bayésienne en R. Il permet de fitter presque n’importe quel type de modèle de régression. Cela comprend tous les modèles que nous avons vu en cours, mais également bien d’autres. Entre autres, des modèles multivariés (i.e., avec plusieurs outcomes), des modèles “distributionnels” (e.g., pour prédire des différence de variance), des generalised additive models, des procesus Gaussiens (Gaussian processes), des modèles issus de la théorie de détection du signal, des mixture models, des modèles de diffusion, des modèles non-linéaires

N’hésitez pas à me contacter pour plus d’informations sur ces modèles ou si vous avez des questions par rapport à vos propres données. Vous pouvez aussi contacter le créateur du paquet brms, très actif en ligne (voir son site). Voir aussi le forum Stan.

References

Efron, B., & Morris, C. (1977). Stein’s paradox in statistics. Scientific American, 236(5), 119–127. https://doi.org/10.1038/scientificamerican0577-119
Gelman, A. (2005). Analysis of variance? Why it is more important than ever. The Annals of Statistics, 33(1), 1–53. https://doi.org/10.1214/009053604000001048
Gelman, A., & Hill, J. (2006). Data analysis using regression and multilevel/hierarchical models. https://doi.org/10.1017/cbo9780511790942
Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao, Y., Kennedy, L., Gabry, J., Bürkner, P.-C., & Modrák, M. (2020). Bayesian workflow. arXiv:2011.01808 [Stat]. http://arxiv.org/abs/2011.01808
Lewandowski, D., Kurowicka, D., & Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of Multivariate Analysis, 100(9), 1989–2001. https://doi.org/10.1016/j.jmva.2009.04.008